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The  quantitative  analysis  and  interpretation  of  button-cell  experiments  usually  depends  upon  assuming 
isothermal  conditions  together  with  uniform  and  known  gas  composition  within  the  gas  compartments. 
An  objective  of  the  present  effort  is  to  develop  computational  tools  to  study  the  validity  of  such  assump¬ 
tions.  A  three-dimensional  computational  fluid  dynamics  (CFD)  model  is  developed  and  applied  to  a 
particular  SOFC  button  cell,  characterizing  the  fluid  flow,  chemistry,  and  thermal  transport.  Results  show 
that  when  inlet  flow  rates  are  sufficiently  high,  button-cell  data  can  be  interpreted  using  the  commonly 
used  assumptions.  However,  when  flow  rates  are  not  sufficient,  the  assumptions  of  uniform  composition 
can  be  significantly  violated.  Additionally,  depending  on  operating  conditions  there  can  be  significant 
temperature  variations  within  the  gas  compartments  and  the  membrane-electrode  assembly. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

The  primary  objective  of  this  paper  is  to  use  a  three-dimensional 
computational  model  to  explore  the  performance  of  a  particular 
solid-oxide  fuel  cell  (SOFC)  button  cell  over  a  range  of  operating 
conditions.  A  supporting  objective  is  to  report  on  the  development 
of  a  three-dimensional  computational  fluid  dynamics  (CFD)  model 
that  extends  the  state  of  the  art,  especially  in  its  ability  to  represent 
elementary  chemistry  and  electrochemistry.  Although  the  study 
considers  a  particular  cell  geometry,  significant  broader  implica¬ 
tions  can  be  inferred  from  the  results. 

Button  cells  play  a  major  role  in  the  development  and  evalu¬ 
ation  of  new  SOFC  materials  and  membrane-electrode  assembly 
(MEA)  architectures.  The  underpinning  models  that  describe  mate¬ 
rial  combinations  and  MEA  architectures  are  usually  developed 
and  validated  with  small  button-cell  experiments.  Quantitative 
predictive  simulations  are  playing  increasingly  important  roles 
in  the  development  of  fuel-cell  stacks  and  systems.  Thus,  it  is 
important  that  laboratory-scale  experimental  observations  can  be 
interpreted  quantitatively  and  unambiguously.  Understanding  the 
influence  of  experimental  configurations  and  operating  conditions 
on  measured  cell  performance  is  essential  to  the  development  and 
validation  of  predictive  models  at  the  system  level. 
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Although  there  are  many  variants  of  button-cell  experiments, 
most  share  some  common  characteristics.  The  button  cell  illus¬ 
trated  in  Fig.  1,  which  represents  an  experiment  developed  by 
Barnett  and  co-workers  [1],  is  typical  and  serves  as  the  focus  of 
the  modeling  in  this  paper.  The  cylindrical  MEA  is  bonded  to  a 
ceramic  support  tube  and  the  entire  assembly  positioned  within 
a  temperature-controlled  furnace.  Fuel  is  supplied  via  a  tube  onto 
the  anode  side  of  the  MEA.  The  cathode  side  is  exposed  to  air, 
which  may  simply  be  the  furnace  air  or  may  be  supplied  via  a  feed 
tube.  The  assembly  shown  in  Fig.  1  is  oriented  horizontally.  Others 
position  the  assembly  vertically.  Button-cell  experiments  should 
usually  be  designed  and  operated  such  that  the  gas-phase  compo¬ 
sition  at  MEA  interfaces  is  known  and  is  spatially  uniform  across 
the  surface  of  the  MEA.  As  a  result,  any  complications  associated 
with  gas-phase  fluid  flow  can  be  removed  from  interpreting  the 
electrochemical  performance.  Depending  upon  gas  flow  rates,  the 
gas-phase  composition  at  the  electrode  faces  can  be  significantly 
different  from  the  composition  entering  through  the  feed  tubes. 

Commercial  CFD  software  can  be  used  directly,  or  modified,  to 
model  aspects  of  fuel-cell  operation,  especially  considering  flow  in 
complex  geometries.  Indeed,  there  is  significant  and  beneficial  use 
of  such  models  to  assist  the  design  and  development  of  cells  and 
stacks.  However,  current  software  and  models  are  usually  limited 
in  their  ability  to  model  details  of  transport,  chemistry,  and  elec¬ 
trochemistry.  For  example,  internal  reforming  chemistry  is  handled 
almost  exclusively  with  two  global  reactions  (steam  reforming  and 
water-gas  shift)  that  use  highly  empirical  rate  expressions.  Charge- 
transfer  chemistry  usually  relies  upon  a  global  Butler-Volmer 
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pre-exponential  factor  for  H2 desorption 
(m2  smol-1) 

specific  surface  area  of  the  active  catalysts  (m-1 ) 
permeability  (m2) 

effective  diffusion  coefficients  (m2  s-1 ) 
binary  diffusion  coefficients  (m2  s-1 ) 
effective  Knudsen  diffusion  coefficients  (m2  s-1 ) 
cell  potential  (V) 

activation  energy  for  H2  desorption  (J  mol-1 ) 
activation  energy  of  electric  conductivity  of  the  elec¬ 
trolyte  (J  mol-1 ) 

activation  energy  for  H2  electrochemical  oxidation 
(Jmol-1) 

activation  energy  for  02  electrochemical  reduction 
(J  mol-1 ) 

reversible  cell  potential  (V) 

Faraday  constant  (Cmol-1 ) 

Grashof  number 

acceleration  of  the  gravity  (m  s-2) 
enthalpy  of  the  gas-phase  mixture  (J  kg-1 ) 
enthalpy  of  the  k  th  species  (J  kg-1 ) 
net  current  density  (Am-2) 
exchange  current  density  (Am-2) 
exchange  current  density  for  H2  oxidation  (Am-2) 
exchange  current  density  for  02  reduction  (Am-2) 
parameter  in  the  expression  of  ig  (Am-2) 
parameter  in  the  expression  of  (Am-2) 
gas-phase  species  mole  flux  (mol  m-2  s-1 ) 
gas-phase  species  mass  flux  (kgm-2  s-1) 
gas-phase  species  mole  flux  within  the  anode 
(molm-2  s-1) 

gas-phase  species  mole  flux  within  the  cathode 
(molm-2  s-1) 

molar  flux  of  H2  (mol  m-2  s-1 ) 
molar  flux  of  H20  (mol  m-2  s-1 ) 
molar  flux  of  02  (mol  m-2  s-1 ) 
number  of  gas-phase  species 
number  of  surface-adsorbed  species 
characteristic  length  (m) 
electrolyte  thickness  (m) 
pressure  (Pa) 
partial  pressure  of  H2  (Pa) 
partial  pressure  of  H20  (Pa) 
partial  pressure  of  02  (Pa) 
parameter  in  the  expression  of  ig  (Pa) 
parameter  in  the  expression  of  ig  (Pa) 
radiative  heat  flux  (J  m-2  s-1 ) 
universal  gas  constant  (J  mol-1  K-1 ) 

Reynolds  number 
pore  radius  (m) 

molar  production  rate  by  surface  reactions 
(molm-2  s-1) 
temperature  (I<) 

MEA  surface  temperature  (K) 
reference  temperature  (I<) 
furnace  temperature  (I<) 
characteristic  velocity  (m  s-1 ) 
gas-phase  velocity  (m  s-1 ) 
molecular  weight  of  H2  (kg  mol-1 ) 


Wk  species  molecular  weight  (kg  mol-1 ) 

W  mean  molecular  weight  (kg  mol-1 ) 

Xk  gas-phase  species  mole  fractions 

Yk  gas-phase  species  mass  fractions 

Greek  letters 

ofa.a  anodic  symmetric  factor  in  the  BV  equation  at  the 
anode 

ac,a  cathodic  symmetric  factor  in  the  BV  equation  at  the 
anode 

aa,c  anodic  symmetric  factor  in  the  BV  equation  at  the 
cathode 

ac<c  cathodic  symmetric  factor  in  the  BV  equation  at  the 
cathode 

volumetric  thermal  expansion  coefficient 
r  surface  site  density  (mol  m-2 ) 

yo  stick  coefficient  of  H2  adsorption 

s  MEA  surface  emissivity 

?7ohm  ohmic  overpotential  (V) 

?7act,a  activation  overpotential  at  the  anode  (V) 

Pact,c  activation  overpotential  at  the  cathode  (V) 

0k  site  fractions  of  surface-adsorbed  species 

A  heat  conductivity  (J  m-1  K  s) 

Aeff  effective  heat  conductivity  (J  m-1  K-1  s-1 ) 

A§  gas-phase  heat  conductivity  (J  m-1 1<-1  s-1 ) 

As  solid-phase  heat  conductivity  (J  m-1  K-1  s-1 ) 
fi  gas-phase  viscosity  (kg  cm-1  s-1 ) 

lik  species  chemical  potentials  (J  mol-1 ) 

/xh2o  chemical  potentials  of  EI20  (J  mol-1 ) 

/xh2  chemical  potentials  of  EI2  (J  mol-1 ) 

/zo2  chemical  potentials  of  02  (J  mol-1 ) 

/x°k  species  standard-state  chemical  potentials  (J  mol-1 ) 

standard-state  chemical  potential  of  EI2  (J  mol-1 ) 
Ph  o  standard-state  chemical  potential  of  EI20  (J  mol-1 ) 

/Xq2  standard-state  chemical  potential  of  02  (J  mol-1 ) 

p  gas-phase  mass  density  (kg  m-3 ) 

g  Stefan-Boltzmann  constant  (W  m-2  K-4) 

r  tortuosity 

0  porosity 


Furnace 


Fig.  1.  Illustration  of  a  nominal  button-cell  system. 
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formulation.  Porous-media  transport  is  handled  in  terms  of  Darcy 
flow.  The  model  reported  herein  relaxes  all  of  these  assumptions, 
and  thus  contributes  to  the  continuing  development  of  high  level 
CFD  software. 

Several  two-  and  three-dimensional  models  have  been  pub¬ 
lished  [2-7]  that  specifically  consider  button-cell  performance. 
Janardhanan  and  Deutschmann  [3]  and  Shi  et  al.  [5,6]  mod¬ 
eled  axisymmetric  (two-dimensional)  button  cells,  but  assumed 
isothermal  conditions.  The  present  three-dimensional,  non- 
isothermal,  model  significantly  extends  these  earlier  contributions. 

2.  Mathematical  models 

The  model  provides  the  capability  for  full  three-dimensional 
fluid  flow,  heat  and  mass  transport,  porous-media  transport  and 
chemistry  throughout  the  gas  compartments  and  the  MEA.  In  the 
open  gas  compartments  (i.e.,  outside  the  porous  electrodes)  the 
flow  is  represented  by  the  Navier-Stokes  equations.  The  chemistry 
and  porous-media  models  are  based  upon  models  developed  by 
Zhu  et  al.  [8,9].  The  solutions  are  constrained  by  boundary  condi¬ 
tions  as  well  as  interface  conditions  between  porous  electrodes  and 
the  gas  compartments. 

2.1.  Gas-phase  transport 

Within  the  gas  compartments,  the  steady-state  mass-,  species-, 
momentum-,  and  energy-conservation  equations  are  stated  gener¬ 
ally  in  vector  form  as 

V-(pV)  =  0,  (1) 

v  •  (pyjfeV)  = -v  -  jfc,  (2) 


/  Kg  \  ^ 

V  .  (phV)  =  V  .  j  A.effVT  -  Y^hkik  +  J2^AsWkhk,  (9) 

\  k=  1  /  k= 1 

sk(0k)  =  0,  (fc  =  1 . JCs)  (10) 

In  these  equations  the  independent  variables  are  the  spatial  coor¬ 
dinates.  The  dependent  variables  are  the  mass  density  p,  velocity  V, 
pressure  p,  mass  fractions  Yk  of  the  I<g  gas-phase  species,  and  cov¬ 
erages  0k  of  the  Ks  surface-adsorbed  species.  Viscous  drag  within 
the  porous  electrodes  is  represented  as  a  Darcy  loss,  where  Bg  is  the 
permeability. 

Within  the  porous  media,  the  gas  and  solid  phases  are  assumed 
to  have  the  same  temperature  T.  The  effective  thermal  conductivity 
Aeff  considers  both  gas-  and  solid-phase  contributions  as 

Xeff  =  #s  +  (l  -0)XS,  (11) 

where  0  is  the  porosity  and  Xs  is  the  solid-phase  thermal  conduc¬ 
tivity.  For  Ni-YSZ,  Xs  «  11  W  m-1  K-1 . 

Heterogeneous  reactions  can  occur  on  the  solid  surfaces  within 
the  porous  anode  structure.  The  molar  production  rates  of  gas- 
phase  species  via  heterogeneous  surface  reaction  are  represented 
as  sk  (Eqs.  (7)  and  (9)).  The  effective  catalytic  specific  surface  area 
(m2  m-3)  is  represented  as  As.  Eq.  (10)  represents  the  fact  that  the 
net  production  rates  ofthe/Cs  surface-adsorbed  species  must  vanish 
at  steady  state.  The  production  rates  of  surface  species  depend  upon 
temperature,  pressure,  gas-phase  composition,  and  the  coverages 
of  surface  species  0k.  Details  about  formulation  of  the  heteroge¬ 
neous  chemistry  may  be  found  in  Kee  et  al.  [10]. 

2.3.  Interface  conditions 


(pW)  =  V  .  rAt(VV  +  (VV)T)-|/i(V-V)l 

-Vp  +  pg 

(3) 

(phV)  =  V  •  (xWT-  V/ifcj*')  , 

(4) 

V  k=  1  J 

Po 

h 

KT^yk/Wk 

(5) 

k=  1 


In  these  equations  the  independent  variables  are  the  spatial  coor¬ 
dinates.  The  dependent  variables  are  the  mass  density  p,  velocity 
V,  pressure  p,  mass  fractions  Yk,  and  temperature  T.  The  nominally 
uniform  pressure  used  in  the  equation  of  state  is  p0.  The  species  dif¬ 
fusive  mass  flux  is  represented  as  jk,  which  is  evaluated  according 
to  ordinary  multicomponent  diffusion.  Thermodynamic  proper¬ 
ties  include  the  mixture  specific  enthalpy  h,  and  species  specific 
enthalpies  hk}  all  of  which  are  known  functions  of  temperature. 
Transport  properties  include  the  gas-phase  dynamic  viscosity  fi 
and  the  thermal  conductivity  X.  Other  variables  are  the  species 
molecular  weights  Wk  and  the  acceleration  of  gravity  g. 

2.2.  Porous-media  transport 

Within  the  porous  electrodes,  the  steady-state  mass-,  species-, 
momentum-,  and  energy-conservation  equations  are  stated  in  vec¬ 
tor  form  as 

V-(pV)  =  0,  (6) 

V  •  (pYkV)  =  -V  ■  jk  +  skAsWk,  (lc=  1 . Kg)  (7) 

V  ■  (pW)  =  -  Vp  +  pg  -  ^V,  (8) 


The  model  does  not  consider  spatial  variations  within  the  thin 
dense  electrolyte  and  the  cathode.  Moreover,  the  anode  charge 
transfer  is  confined  to  a  thin  layer  at  the  interface  between  the 
dense  electrolyte  and  the  porous  anode  structure.  In  other  words, 
all  the  electrochemistry  and  ion  transport  is  confined  to  an  infinitely 
thin  layer,  and  thus  the  behavior  is  represented  in  terms  of  inter¬ 
face  conditions.  This  approximation  is  based  upon  the  results  of 
previously  published  models  that  spatially  resolve  the  thin  charge- 
transfer  region,  showing  that  charge  transfer  is  accomplished 
within  a  few  tens  of  microns  [9,11  ]. 

Consider  first  the  interface  between  the  fuel  compartment  and 
the  porous  anode.  The  solutions  to  Eqs.  (l)-(9)  are  constrained  to 
require  that  the  fluxes  of  mass,  species,  momentum,  and  energy 
are  continuous  across  the  interface.  Additionally,  radiation  heat 
transfer  from  the  face  of  the  porous  electrode  is  found  to  have  a 
significant  impact  on  the  MEA  temperature. 

Anodic  charge-transfer  chemistry  proceeds  at  the  interface 
between  the  dense  electrolyte  and  the  porous  anode  structure.  The 
result  is  a  net  molar  flux  of  gas-phase  species  into  and  out  of  the 
anode  pore  spaces  as  represented  as  Jk.  Cathodic  charge-transfer 
chemistry  proceeds  at  the  interface  between  the  thin  cathode  and 
the  air  compartment.  In  this  case,  the  net  flux  of  oxygen  between 
the  air  and  the  cathode  is  represented  as  Jk. 

In  addition  to  the  species  fluxes,  the  model  must  also  represent 
other  important  processes  within  the  thin  membrane-electrode 
assembly.  Irreversible  polarization  losses  associated  with  the  elec¬ 
trochemical  charge-transfer  chemistry  manifest  themselves  as 
heat.  The  heating  rate  is  evaluated  from  the  enthalpy  fluxes  to  and 
from  the  thin  MEA  layer.  The  electrical  power  density  produced 
by  the  cell  is  represented  as  i£cell,  where  i  is  the  electrical  cur¬ 
rent  density  and  Ecell  is  the  cell  operating  potential.  This  electrical 
energy  is  subtracted  from  the  thermal  energy  balance.  There  is  also 
a  thermal-radiation  loss  from  the  cathode  side  of  the  MEA. 
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2.4.  Boundary  conditions 

Boundary  conditions  for  the  gas  flows  are  relatively  straightfor¬ 
ward.  The  composition,  velocity,  and  temperature  of  the  fuel  inlet  is 
specified.  The  fuel-side  outflow  boundary  is  specified  by  specifying 
the  exit  pressure.  The  particular  reactor  used  in  this  study  uses  an 
unforced  air  supply  on  the  cathode  side.  Thus  the  boundary  condi¬ 
tion  is  specified  as  the  composition  of  air  at  atmospheric  pressure 
and  at  furnace  temperature. 

There  are  several  solid  surfaces  in  the  model,  not  including  the 
porous  media.  These  include  the  fuel  feed  tube  and  the  cylindrical 
structure  that  supports  the  button  cell.  All  these  solid  surfaces  are 
taken  to  be  isothermal  at  a  specified  furnace  temperature  and  the 
velocity  vanishes  on  all  the  solid  surfaces.  Additionally,  the  solid 
surfaces  are  assumed  to  be  chemically  inert,  such  that  the  species 
fluxes  normal  to  the  surface  vanish. 

2.5.  Thermal  radiation 

Despite  the  fact  that  temperatures  and  temperature  differences 
are  relatively  low,  there  is  the  possibility  of  thermal  radiation  from 
the  surfaces  of  the  MEA.  Although  the  thermal  radiation  is  weak,  the 
low  flow  rates  also  cause  a  weak  convective  interaction  between  the 
gas  compartments  and  the  MEA.  The  present  study  uses  a  relatively 
simple  radiation  model  in  which  the  surfaces  of  the  MEA  radiate  to 
a  black-body  enclosure  at  the  nominal  furnace  temperature.  That 
is,  the  radiative  flux  from  the  MEA  surfaces  is  evaluated  as 

<'ad=^(Cf-ri)>  (12) 

where  a  is  the  Stefan-Boltzmann  constant,  Tsurf  is  the  MEA  sur¬ 
face  temperature,  and  is  the  furnace  temperature.  In  the  present 
study,  the  emissivity  of  the  MEA  surfaces  is  taken  to  be  €  =  0.8.  This 
value  is  larger  than  the  emissivities  of  smooth  Nickel  or  Ziconia 
since  the  anode  pore  spaces  trap  incident  radiation. 

2.6.  Diffusion  model 

The  diffusion  mass-flux  vector,  jk,  is  modeled  in  both  gas-phase 
fuel  and  air  chambers,  as  well  as  the  porous  anode,  using  a  Fickian 
model: 


The  results  reported  in  this  paper  assume  that  the  anode  porous 
matrix  is  homogenous  and  isotropic.  However  the  model  itself  is 
capable  of  spatially  varying  and  anisotropic  properties,  which  may 
be  needed  to  represent  certain  functionally  graded  electrode  struc¬ 
tures. 

2.7.  Heterogeneous  chemistry 

Heterogeneous  chemistry  (e.g.,  reforming,  partial  oxidation, 
etc.)  occurs  on  active  catalyst  surfaces  within  the  porous  anode 
structure.  Although  the  models  and  software  are  quite  general  and 
accept  alternative  reaction  mechanisms,  the  mechanism  used  in 
the  present  study  considers  reaction  within  Ni-YSZ  cermet  struc¬ 
tures.  This  reaction  mechanism,  which  has  been  validated  and 
documented  previously  [3,8],  consists  of  42  reactions  among  6 
gas-phase  species  and  12  surface-adsorbed  species.1  The  surface 
reaction  rates  have  Arrhenius  temperature  dependence,  and  some 
reactions  have  activation  energies  that  are  functions  of  the  surface 
coverages. 

This  catalytic  reaction  mechanism  spans  the  conditions  between 
steam  reforming,  autothermal  reforming,  dry  reforming  (i.e., 
reforming  of  CH4  by  CO2),  and  catalytic  partial  oxidation.  This 
approach  is  quite  different  from,  and  more  capable  than,  most 
published  SOFC  models,  which  usually  incorporate  two  global 
reactions  for  steam  reforming  and  water  gas  shift.  The  global 
approach  requires  empirical  rate  expression  for  the  global  reac¬ 
tions.  The  present  elementary  mechanism  does  not  directly  include 
a  steam-reforming  or  a  water-gas-shift  reaction.  Rather,  the  pro¬ 
cesses  are  represented  by  combinations  of  elementary  reactions. 
This  fundamental  approach  enables  predictive  capabilities  over 
wide  operating  conditions,  spanning  steam  reforming,  autother¬ 
mal  reforming,  and  partial  oxidation.  Water-gas-shift  chemistry  is 
also  fully  represented  within  the  reaction  mechanism. 

2.8.  Electrochemistry 

The  operating  cell  potential  £cell  can  be  written  in  terms  of  over¬ 
potentials  as 

ECell  —  £rev  —  ^?act,a(0  —  l^act,c(0l  —  ^ohm(0-  (17) 


jk  =  ~pDekVYk.  (13) 

For  the  fuel  and  air  chambers,  the  effective  diffusion  coefficient 
is  the  mean  bulk  diffusion  coefficient,  D™,  calculated  as, 


D ™ 


l-Xfc 

Kg 


(14) 


where  Xk  are  the  mole  fractions,  Vkj  are  the  binary  diffusion  coef¬ 
ficients  determined  from  kinetic  theory  [10],  and  y  =  1. 

Within  the  porous  anode,  the  effective  mixture-averaged  diffu¬ 
sion  coefficients  are  evaluated  using  the  Bosanquet  formula, 


1  _  1  1 
DTW  +  Df 


(15) 


The  mean  bulk  diffusion  coefficients  D™  are  evaluated  from  Eq.  (14) 
with  y  set  to  0/r,  where  0  is  the  porosity  and  r  is  the  tortuosity. 
The  Knudsen  diffusion  coefficient  in  Eq.  (15)  is  calculated  as 
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In  Eq.  (17),  Erev  is  the  reversible  (Nernst)  potential  across  the 
dense  electrolyte,  ijohm  is  the  ohmic  overpotential,  and  ?7act,a  and 
?7act,c  are  the  activation  overpotentials  at  the  anode-electrolyte  and 
cathode-electrolyte  interfaces,  respectively.  All  the  overpotentials 
depend  upon  the  current  density.  Because  species  transport  and 
chemistry  are  fully  resolved  within  the  anode,  a  model  for  concen¬ 
tration  overpotential  is  not  needed. 

In  a  typical  porous,  cermet,  composite,  electrode,  electrochem¬ 
ical  reactions  occur  within  a  very  thin  region  (tens  of  microns) 
adjacent  to  the  dense  electrolyte  [9,11  ].  With  a  button  cell  that  is  on 
the  order  of  1  mm  thick,  it  is  not  necessary  to  resolve  the  charge- 
transfer  region  spatially.  Rather,  it  is  sufficient  to  approximate  the 
charge-transfer  behavior  as  though  it  proceeds  at  a  surface.  In 
this  approach,  geometric  factors  such  as  the  three-phase-boundary 
lengths  are  subsumed  into  expression  for  the  exchange  current 
densities  [8]. 

The  results  in  this  paper  are  based  upon  assuming  that  electro¬ 
oxidation  of  the  fuel  proceeds  entirely  via  H2.  It  is  known  that  CO, 
and  possibly  hydrocarbons,  can  participate  in  charge-transfer  reac¬ 
tions  [12,13].  However,  the  charge-transfer  rates  are  greatly  higher 
for  H2.  Moreover,  catalytic  water-gas-shift  processes  (globally, 


where  rp  is  the  mean  pore  radius. 


1  The  mechanism  is  available  for  download  at  http://www.detchem.com. 
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CO  +  H20  ^  C02  +  H2)  readily  convert  CO  to  H2.  Consequently  it 
is  a  reasonable  to  use  the  widely  used  approximation  that  charge 
transfer  proceeds  through  H2.  As  ongoing  research  reveals  the 
information  needed  to  establish  the  fundamental  reaction  path¬ 
ways  and  rates  for  parallel  charge-transfer  processes,  they  can  be 
incorporated  into  models  such  as  the  ones  discussed  here.  However, 
assuming  H2charge  transfer  is  entirely  sufficient  for  the  purposes  of 
the  present  paper,  which  focuses  upon  thermal  and  flow  behaviors. 

Assuming  that  electro-oxidation  of  the  fuel  proceeds  domi¬ 
nantly  via  H2,  the  reversible  potential  can  be  calculated  from  the 
Nernst  equation  as 


E,pv  -  A  ( 
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RT  , 


A*h2  +  2> ^>2  -  **h2o)  +  Jf  ln 


Wapff; 

PH20,a 


(18) 


where  F  is  the  Faraday  constant,  /x£  are  standard-state  chemical 
potentials,  pk  are  species  partial  pressures  (measured  in  atmo¬ 
spheres).  Subscripts  a  and  c  denote  anode-side  and  cathode-side, 
respectively.  Any  small  losses  due  to  electrical  current  leakage 
through  the  dense  electrolyte  are  neglected. 

Zhu  et  al.  [8]  derived  and  validated  a  modified  Butler-Volmer 
formulation,  which  is  based  upon  a  multi-step  H2  electrochem¬ 
ical  oxidation  process  with  a  single  rate-limiting  step.  At  the 
anode-electrolyte  interface,  the  current  density  is  represented  as 
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where  zg  is  the  anode  exchange  current  density,  and  aa,a  and  ac,a 
are  the  anodic  and  cathodic  symmetry  parameters,  respectively. 
The  anode  exchange  current  density  is  calculated  as 
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where  £h2  is  the  apparent  activation  energy  and  i*  is  an  empir¬ 
ical  constant  that  is  adjusted  to  fit  experiments  at  a  reference 
temperature  Tref.  The  term  p *  is  derived  for  a  steady  hydrogen 
adsorption-desorption  balance: 


Ph2  = 


Ad(,sr2  ^2  jtRTWh, 
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exp 


(21) 


where  the  pre-exponential  factor  Ades,  the  activation  energy  Edes, 
the  surface  site  density  r,  and  the  sticking  coefficient  y0  are  con¬ 
stants  and  Wh2  is  the  molecular  weight  of  atomic  hydrogen. 

At  the  cathode-electrolyte  interface,  the  modified 
Butler-Volmer  equation  is 


=  exp  (^g^)  -  exp  (-^gactc)  . 
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The  cathode  exchange  current  density  is  evaluated  using 
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where  Eq2  is  the  apparent  activation  energy  for  the  electrochemical 
reduction  of  02,  z*2  is  an  adjustable  constant,  and  p*  2  is  written  in 
Arrhenius  form.  Further  details  are  found  in  Zhu  et  al.  [8]. 

The  ohmic  overpotential  (Eq.  (17))  is  modeled  in  terms  of  02-ion 
transport  resistance  within  the  dense  YSZ  electrolyte, 


P°hm  CTo7'-leXp(-Eel/R7')'  ( 

In  this  expression  Lel  is  the  electrolyte  thickness  and  the  denomi¬ 
nator  is  the  temperature-dependent  ion  conductivity. 


The  current  density  i  is  determined  by  substituting  Eqs. 
( 18)-(24)  into  Eq.  (17).  The  species  molar  fluxes  due  to  electrochem¬ 
istry,  J£l  are  determined  from  the  current  density  as  J®1  a  =  -i/2F, 

iH20,a  =  '72fandio2,c  =  -'74F. 


3.  Numerical  method 

The  conservation  equations  are  discretized  using  an  unstruc¬ 
tured  finite-volume  method  with  second-order  reconstruction 
of  the  convective  and  diffusive  fluxes  for  all  equations  [14]. 
This  approach  enables  the  incorporation  of  three-dimensional 
geometric  complexity,  using  geometry-definition  and  automated 
mesh-discretization  software  (e.g.,  Gambit). 

At  all  materials  interfaces,  there  are  finite-volume  cells  that  abut 
the  interface.  In  a  finite-volume  setting,  the  conservation  equa¬ 
tions  are  represented  in  an  integral  sense  for  each  finite  volume. 
Away  from  material  interfaces,  fluxes  leaving  one  cell  must  be 
identical  to  the  fluxes  entering  the  neighboring  cell.  At  material 
interfaces,  however,  certain  physical  or  chemical  behaviors  can 
manifest  themselves  as  sources  of  mass,  energy,  etc.  For  exam¬ 
ple,  thermal  radiation  at  the  interface  between  a  porous  material 
and  the  neighboring  gas  flow  represents  a  loss  (or  gain)  in  ther¬ 
mal  energy  for  the  finite  volumes  that  span  the  interface.  A  thin 
(i.e.,  not  spatially  resolved)  MEA  is  a  more  complex  interface.  The 
charge-transfer  chemistry  causes  gas-phase  species  fluxes  into  and 
out  of  neighboring  finite  volumes.  Moreover,  there  are  important 
thermal  consequences.  Polarization  losses  associated  with  charge- 
transfer  chemistry  and  Joule  heating  associated  with  ion  current 
through  the  dense  electrolyte  are  sources  of  thermal  energy  at  the 
interface.  However,  net  electrical  power  delivered  to  an  external 
load  (i.e.,  zECeii)  does  not  contribute  to  the  thermal  balance. 

The  system  of  equations  is  solved  using  a  segregated 
SiMPLEalgorithm,  in  which  the  mass-continuity  equation  (Eqs.  (1) 
and  (6))  are  recast  as  a  pressure-correction  equations  and  solved 
together  with  the  momentum  equations  (i.e.,  Eqs.  (3)  and  (8))  [15]. 
The  species-continuity  and  thermal-energy  equations  are  similarly 
linearized.  The  inner  linear  problem  is  solved  using  an  algebraic 
multigrid  (AMG)  method  [16].  An  outer  iteration  is  used  to  solve 
the  full  nonlinear  problem. 

The  software  is  written  in  C++  using  object-oriented  pro¬ 
gramming  techniques,  and  interfaced  with  the  scriptable 
PYTHONlanguage.  The  Cantera2  software,  which  is  also  written  in 
C++,  is  used  to  evaluate  thermodynamic  quantities  (enthalpies, 
heat  capacities,  Nernst  potentials),  transport  properties  (viscosi¬ 
ties,  thermal  conductivities,  species  diffusivities)  and  reaction 
rates  [17]. 

3.1.  Verification 

The  software  that  implements  the  model  is  large  and  complex, 
providing  three-dimensional  geometrical  generality  and  the  incor¬ 
poration  of  extensive  chemical  and  physical  processes.  Thus,  as 
in  all  major  software  development,  comprehensive  verification  is 
an  important  aspect  of  effort.  In  this  context,  verification  is  con¬ 
cerned  primarily  with  assuring  that  the  software  solves  accurately 
the  underlying  physical  model.  Software  verification,  per  se,  is  less 
concerned  with  validating  the  physical  and  chemical  formulations, 
which  usually  rely  upon  experimental  measurements  and  their 
interpretation.  Of  course,  validating  the  underpinning  physics  and 
chemistry  is  also  an  important  objective. 


2  See  http://www.cantera.org. 
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Fig.  2.  Gas-phase  mole  fractions  along  the  centerline  of  an  SOFC  planar  fuel  channel. 
Predictions  from  the  current  CFD  model  are  shown  as  solid  lines  and  predictions 
from  a  plug-flow  model  [8]  are  shown  as  dashed  lines. 

The  object-oriented  program  structure  enables  the  independent 
testing  and  validation  of  all  the  sub-models  (thermodynamic  and 
transport  property  evaluations,  chemical  kinetics,  etc.).  Addition¬ 
ally,  the  present  model  has  been  tested  against  other  independently 
developed  models  that  solve  similar  problems.  Such  model-to- 
model  comparisons  are  an  important  aspect  of  model  development 
and  validation.  For  example,  the  three-dimensional  model  was 
tested  by  comparing  a  planar  SOFC  channel  simulation  with  results 
published  by  Zhu  et  al.  [8],  using  conceptually  identical,  but  inde¬ 
pendently  programmed,  submodels  and  model  parameters. 

Fig.  2  compares  gas-phase  mole  fractions  along  the  anode- 
channel  centerline.  The  present  CFD  model  fully  resolves  the 
channel  fluid  flow  while  the  model  by  Zhu  et  al.  [8]  assumes  plug 
flow.  The  CFD  model  enables  two-dimensional  species  transport 
within  the  anode  while  Zhu  et  al.  neglect  axial  diffusive  trans¬ 
port.  Despite  these  differences  in  the  underlying  models,  the  results 
agree  well.  These  comparisons  serve  to  validate  the  current  CFD 
model  (especially  the  chemistry  and  electrochemistry  submodels) 
as  well  as  the  low-dimensional  approximations  in  the  model  by  Zhu 
et  al.  It  should  also  be  noted  that  chemistry  aspects  of  the  Zhu  et  al. 
model  were  developed  and  validated  using  experimental  data. 

4.  Specific  button  cell 

The  button-cell  geometry  used  in  the  present  effort  represents 
an  experimental  facility  developed  by  Barnett  et  al.  [1].  As  illus¬ 
trated  in  Fig.  1 ,  the  cell  consists  of  an  inlet  fuel  tube  (6.5  mm  outside 
diameter)  with  four  cylindrical  fuel  channels  (1.8  mm  diameter). 
The  fuel  tube  is  positioned  coaxially  with  an  alumina  tube  (9.5  mm 
inside  diameter)  that  supports  the  MEA.  The  end  of  the  fuel  tube 
is  positioned  3  mm  from  the  MEA.  Anode  exhaust  gases  flow  out 
axially  through  the  annular  space  between  the  fuel  tube  and  the 
MEA  support  tube. 

The  cathode  side  of  the  MEA  is  exposed  to  air  within  the  furnace. 
As  illustrated  in  Fig.  1  there  is  no  forced-air  supply.  The  assembly  is 
mounted  horizontally  in  a  furnace  whose  walls  are  maintained  at 
uniform  temperature.  Thus,  three-dimensional  effects  associated 
with  natural  convection  are  a  potential  concern. 

The  MEA  is  fabricated  using  a  Ni-YSZ  porous  anode  that  has 
a  thickness  of  0.7  mm  and  a  diameter  of  15.2  mm.  The  dense 
electrolyte  is  10  |xm  of  YSZ,  and  the  cathode  is  40  [xm  of  porous 
YSZ-LSM.  The  cathode  is  only  8  mm  in  diameter,  so  does  not  span 
the  entire  diameter  of  the  anode  disk.  Thus,  only  the  inner  section 
of  the  MEA  is  electrochemically  active. 

The  geometry  is  symmetric  about  a  90°  sector.  Flowever,  because 
the  cell  is  mounted  horizontally,  a  180°  sector  is  modeled  to  include 


Fig.  3.  Button-cell  geometry  and  mesh,  showing  fuel-tube  inlet  and  outlet,  walls, 
MEA  and  air  plenum. 

natural  convection.  Fig.  3  shows  some  of  the  unstructured  mesh 
faces.  In  total  the  mesh  is  comprised  of  77,336  cells. 

The  four-channel  fuel  tube  is  somewhat  unusual;  most  button¬ 
cell  experiments  use  a  single  inlet  tube.  While  the  present  model 
uses  the  exact  geometry  for  a  particular  experiment  [1  ],  the  results 
depend  only  weakly  on  the  inlet  tube  details.  Flowever,  the  results 
depend  greatly  upon  inlet  flow  rates.  It  may  be  noted  that  the 
present  study  uses  a  relatively  small  button  cell.  Nevertheless,  the 
general  observations  and  insights  from  the  study  are  applicable  to 
larger  cells. 

5.  Model  parameters 

Models  such  as  the  one  here  necessarily  involve  a  number 
of  parameters.  Some,  such  as  physical  dimensions,  are  measured 
directly.  Flowever,  as  a  practical  matter,  there  are  model  param¬ 
eters  that  must  be  established  empirically.  These  usually  involve 
electrochemical  reaction  rates  and  available  three-phase  bound¬ 
ary  area,  which  are  not  measured  independently.  The  parameters 
shown  in  Table  1  are  determined  to  fit  experimental  data  (Fig.  4) 
using  humidified  hydrogen  fuel  at  nominal  furnace  temperatures  of 
700, 750  and  800  °C.  The  fuel  flow  rates  are  set  to  be  sufficiently  high 
that  depletion  and  dilution  are  negligible  (200  seem).  By  flooding 
the  cell,  the  charge-transfer  parameters  can  be  determined  without 
concern  for  fluid-mechanical  effects  in  the  anode  chamber.  Once 
established,  the  parameters  are  held  fixed  for  all  other  simulations 
in  this  study.  Note  that  at  200  seem,  for  the  highest  current  density 
simulated  (0.1V),  the  computed  fuel  depletion  is  45%.  In  contrast, 
at  the  lowest  flow  rate  simulated,  10  seem,  the  fuel  depletion  is 
73%. 

Fig.  4  a  compares  measured  and  predicted  polarization  charac¬ 
teristics  using  the  full  three-dimensional  model.  It  is  evident  that 
the  set  of  model  parameters  produces  an  excellent  representation 
of  the  data.  Important  fitted  parameters  include  the  exchange  cur¬ 
rent  factors  and  z^,  and  the  activation  energies  Eh2  and  Eo2- 
Although  used  as  fitting  parameters  to  match  data  in  Fig.  4,  all 
the  parameters  are  physically  reasonable.  Nevertheless,  it  must 
be  acknowledged  that  the  set  of  parameters  may  not  be  entirely 
unique. 
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Table  1 

Physical  and  chemical  parameters  for  the  MEA  structure. 


Parameters 

Value 

Units 

Anode 

Thickness  (La) 

700 

fxm 

Porosity  (0) 

0.4 

Tortuosity  (r) 

4.8 

Mean  pore  radius  (rp) 

0.125 

pan 

Particle  diameter  (dp) 

1.0 

|jim 

Specific  catalyst  area  (As) 

10000 

cm-1 

Apparent  activation  energy  (Eh2  ) 

110.0 

kj  mol-1 

Reference  temperature  (Tref) 

1023.0 

K 

Exchange  current  factor  (i^ ) 

1.25 

Acm~2 

Anodic  symmetry  factor  (aa>a) 

1.5 

Cathodic  symmetry  factor  (ac,a) 

0.5 

pf,  pre-exponential  (Ades ) 

5.59E19 

sem2  mol-1 

p* j  activation  Energy  (Ede s) 

88.12 

kj  mol-1 

p* j2  surface  site  density  (T) 

2.6E-9 

mol  cm-2 

p* j2  sticking  probability  (y0) 

0.01 

Cathode 

Apparent  activation  energy  ( Eq2  ) 

110.0 

kj  mol-1 

Reference  temperature  (Tref) 

1023.0 

I< 

Exchange  current  factor  (i*  ^ ) 

1.25 

Acm~2 

Anodic  symmetry  factor  (aa>c) 

0.875 

Cathodic  symmetry  factor  (ac,c) 

0.125 

Electrolyte:  crel  =  cr0/T  exp(-Ee\/RT ) 

Thickness  (Le i) 

10 

pan 

Activation  energy  (Eel) 

80.11 

kj  mol-1 

Ion  conductivity  pre-factor  (cr0) 

1.786E5 

serrr1 

i - 1 - 1 - 1 - r 


0  12  3 


Current  density  (A/cm2) 

Fig.  4.  Experimental  and  simulated  polarization  curves  for  a  button  cell  operating 
with  flooded  (200  seem)  humidified  H2.  The  upper  panel  shows  predictions  with  a 
full  model  that  includes  the  energy  equations  and  thus  temperature  distributions. 
Model  parameters  were  established  with  the  full  model.  The  lower  panel  shows 
predicted  polarization  when  assuming  isothermal  conditions  at  the  nominal  furnace 
temperature. 


Button-cell  data  are  generally  reported  at  certain  temperatures, 
usually  inferring  isothermal  operation.  However,  as  discussed  sub¬ 
sequently,  temperature  can  vary  significantly  throughout  the  MEA, 
especially  at  high  current  density.  The  MEA  temperatures  are  usu¬ 
ally  greater  than  the  furnace-wall  temperature  due  to  heat  released 
by  electrochemical  reactions  and  Ohmic  heating.  Fig.  4  b  shows  pre¬ 
dicted  polarization  behavior  using  the  three-dimensional  model, 
but  without  solving  the  energy  equations  and  maintaining  isother¬ 
mal  conditions  of  700,  750  and  800  °C.  By  comparing  predictions 
from  the  full  model  with  those  under  isothermal  conditions,  it 
is  evident  that  temperature  variations  within  the  MEA  can  have 
a  significant  effect  on  polarization  behavior.  These  results  show 
that  assuming  isothermal  conditions  when  interpreting  button-cell 
experiments  can  be  misleading. 

6.  Performance  on  hydrogen 

Fig.  5  shows  experimentally  measured  characteristics  for  the 
button  cell  operating  on  humidified  H2,  but  at  different  fuel  flow 
rates.  At  flow  rates  above  200  seem,  the  fuel  depletion  is  found 
to  be  negligible  and  the  performance  is  essentially  independent  of 
flow  rate.  At  low  flow  rates  the  effects  of  fuel  depletion  and  dilu¬ 
tion  can  be  significant,  especially  at  high  current  density.  However, 
the  polarization  behavior  is  surprisingly  weakly  affected  by  the 
flow  rate.  Fig.  5  indicates  the  low  flow  rate  (25  seem)  causes  lower 
voltages,  even  at  open  circuit.  However,  the  open-circuit  potential 
should  be  independent  of  flow  rate.  The  approximately  30  mV  off¬ 
set  for  the  25  seem  data  is  attributable  to  minor  gas  leaks  in  the 
experimental  setup.  The  simulations  neglect  this  leakage  effect  at 
all  flow  rates.  With  inlet  fuel  flow  rates  of  25  seem,  gas  chromatog¬ 
raphy  suggests  that  approximately  1%  oxygen  is  getting  into  the  fuel 
flow.  This  effect  is  too  small  to  have  any  impact  on  cell  performance 
at  high  flow  rates. 

Fig.  6  shows  model  predictions  for  low  (25  seem)  and  high 
(200  seem)  flow  rates.  In  both  cases  the  cell  operating  with  a 
fuel  mixture  of  97%  H2,  3%  H20.  The  wall  temperatures  and  flow 
inlet  temperatures  are  fixed  at  800  °C.  The  operating  potential  is 
FCeii  =  0.5  V.  To  help  visualize  the  flow,  path  lines  originating  at  the 
inlets  are  superimposed  on  the  mole-fraction  and  temperature  con¬ 
tours  in  Fig.  6.  Because  the  flow  field  is  essentially  axisymmetric 
under  these  operating  conditions,  contours  on  half  the  symme¬ 
try  plane  are  shown.  Nevertheless,  the  simulations  themselves  are 
three-dimensional. 

At  low  fuel  flow  rates,  it  is  clear  that  the  gas  composition  within 
the  anode  compartment  varies  significantly  from  the  inlet  compo¬ 
sition.  At  sufficiently  high  flow  rates,  however,  the  composition  at 


Current  density  (A/cm2) 


Fig.  5.  Measured  polarization  behavior  at  different  fuel  flow  rates.  In  all  cases,  the 
fuel  is  humidified  H2  and  the  nominal  temperature  is  800  °C. 


130 


G.M.  Goldin  et  al.  /  Journal  of  Power  Sources  187  (2009)  123-135 


Fig.  6.  Contour  plots  illustrating  model  solutions  for  the  cell  operating  with  a  fuel 
mixture  of  97%  H2,  3%  H20,  operating  potential  of  0.5  V,  and  wall  temperature  of 
800  °C.  Solutions  are  shown  for  fuel  flow  rates  of  25  and  200  seem.  The  composition 
contours  are  mole  fractions.  Path  lines  are  superimposed  on  the  contours  to  show 
flow  patterns. 


the  interface  between  the  anode  gas  compartment  and  the  face  of 
the  MEA  is  close  to  that  of  the  inlet  flow.  Interpreting  experimental 
results  usually  depends  upon  knowing  the  gas-phase  composition 
within  the  anode  compartment.  Thus,  it  is  desirable  to  maintain 
sufficiently  high  flow  rates  such  that  there  is  little  fuel  depletion 
and  the  composition  is  nearly  that  of  the  feed  gases. 

Fig.  7  shows  gas-phase  H2  mole-fraction  profiles  along  the 
button-cell  axis,  both  in  the  anode  compartment  and  within  the 
anode  pore  space.  Solutions  are  shown  for  flow  rates  of  200  and 
25  seem.  The  axial  position  of  x  =  0  is  at  the  interface  between 
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Fig.  7.  H2  mole  fraction  along  the  button-cell  axis  in  the  fuel  chamber  and  within 
the  anode  pore  space.  Solutions  are  shown  for  two  fuel  flow  rates:  (a)  200  seem  and 
(b)  25  seem. 


the  porous  anode  and  the  gas  in  the  anode  chamber.  The  right¬ 
most  position  at  z  =  0.7  mm  is  at  the  anode-electrolyte  interface. 
The  cell  is  operating  with  humidified  hydrogen  (i.e.,  97%  H2  and 
3%  H20).  It  is  evident  that  there  is  considerable  hydrogen  varia¬ 
tion,  especially  within  the  porous  anode.  The  variation  is  greatest 
at  low  operating  potential,  where  hydrogen  is  being  consumed  at 
higher  rates  to  support  the  electrochemical  charge  transfer.  The  sig¬ 
nificantly  higher  gradients  within  the  anode  pore  space  are  due  to 
comparatively  high  transport  resistance  in  the  porous  media. 

At  flow  rates  of  200  seem  (Fig.  7a),  there  is  relatively  little  fuel 
depletion  in  the  anode  compartment.  The  H2mole  fraction  remains 
greater  than  90%,  even  at  the  high  current  densities  associated 
with  low  operating  potential.  Under  these  conditions,  the  hydro¬ 
gen  profiles  are  not  affected  greatly  by  the  isothermal  assumption. 
Although  there  is  no  general  rule  for  determining  an  acceptable 
level  of  fuel  depletion,  maintaining  fuel  composition  within  a  few 
percent  of  the  inlet  composition  is  certainly  reasonable.  As  illus¬ 
trated  by  Fig.  7  a,  inlet  flow  of  200  seem  in  this  reactor  appears 
to  be  sufficiently  high  even  at  very  low  operating  potentials.  This 
trend  is  confirmed  by  the  asymptotic  behavior  of  the  experimental 
polarization  curves  at  high  flow  rates  (Fig.  5). 

At  a  relatively  low  fuel  flow  rate  of  25  seem  (Fig.  7b),  there  is 
significant  fuel  depletion.  Especially  at  low  operating  potentials, 
the  H2 composition  in  the  anode  chamber  is  significantly  below  the 
specified  inlet  composition  of  97%.  As  can  also  be  seen  in  Fig.  6, 
at  low  flow  rates  there  is  very  significant  upstream  diffusion  well 
into  the  feed  tube.  Quantitative  interpretation  of  button-cell  exper¬ 
iments  at  these  low  flow  rates  would  depend  upon  a  quantitative 
understanding  of  fuel  depletion. 

In  addition  to  depletion  along  the  centerline  axis,  there  are  also 
radial  composition  variations.  Fig.  8  shows  the  radially  averaged  H2 
mole  fraction  at  the  interface  between  the  anode  and  fuel  chamber 
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Fig.  8.  Average  H2  mole  fraction  at  the  interface  between  the  anode  and  fuel  cham¬ 
ber  as  functions  of  fuel  flow  rate  and  cell  potential.  In  all  cases,  the  inlet  fuel  is  97% 
H2. 

as  functions  of  inlet  fuel  flow  rate  at  various  cell  potentials.  As 
current  density  increases  (decreasing  cell  potential),  more  fuel 
is  utilized  and  the  H2  depletion  increases.  Similarly,  as  the  inlet 
fuel  flow  rate  decreases,  there  is  less  H2  to  mix  with  the  reaction 
product  H20,  and  H2  dilution  increases.  For  the  most  extreme  case 
shown  in  Fig.  8(£cell  =  0.1  V  and  10  seem),  the  average  FI2  mole 
fraction  at  the  anode  interface  decreases  to  around  60%,  which  is 
substantially  smaller  than  the  inlet  mole  fraction  of  97%. 

Fig.  6  c  shows  that  the  maximum  temperatures  within  the  MEA 
are  some  10°  C  higher  than  the  wall  temperatures  and  the  inlet 
gas  temperatures.  Fig.  9  shows  temperature  profiles  through  the 
gas  compartments  and  the  MEA  structure.  The  highest  tempera¬ 
tures  are  always  near  dense  electrolyte,  where  the  charge-transfer 
chemistry  is  active  and  polarization  losses  manifest  themselves  as 
heat  sources.  At  low  operating  potential  (high  current  density),  the 
MEA  temperatures  are  nearly  25 0  C  higher  than  the  furnace  wall 
temperature.  The  temperature  profiles  are  found  to  depend  only 
weakly  upon  the  fuel  rates. 

It  is  interesting  to  note  that  if  radiation  from  the  MEA  sur¬ 
face  is  neglected,  the  predicted  MEA  temperatures  are  significantly 
greater.  Without  radiation,  the  predicted  peak  MEA  temperature  is 
nearly  820°  Cat  Ecell  =  0.5  V,  and  nearly  840 0  C  at  £cell  =  0.1  V.  It  is 
also  interesting  to  note  that  the  temperature  profiles  in  the  gas  com¬ 
partments  are  nearly  linear.  Such  profiles  show  that  heat  transfer 
through  the  gas  phase  is  dominated  by  conduction,  not  convection. 
This  explains  why  the  temperature  profiles  depend  weakly  upon 
flow  fuel  rates. 
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Fig.  9.  Centerline  temperature  profiles  through  the  gas  compartments  and  the  MEA 
structure.  In  all  cases,  the  cell  is  operating  on  an  inlet  flow  rate  of  200  seem  of 
humidified  H2.  The  wall  temperatures  and  inlet  gases  are  at  800  °C. 
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Fig.  10.  Predicted  polarization  characteristics  for  the  button  cell  operating  at  800° 
C  with  a  fuel  mixture  of  97%  H2,  3%  H20,  but  with  fuel  flow  rates  of  either  25  or 
200  seem. 


It  is  clear  that  low  flow  rates  have  a  significant  effect  on  fuel 
composition.  Interestingly,  however,  the  flow-rate  effects  on  polar¬ 
ization  appear  to  be  weak  (Fig.  5).  Fig.  10  shows  model-predicted 
polarization  behavior  for  25  and  200  seem  fuel  flow  rates  (corre¬ 
sponding  to  the  calculations  shown  in  Fig.  6).  Despite  the  significant 
variations  in  fuel  depletion,  the  predicted  effect  on  polarization  is 
very  small.  This  weak  influence  is  partially  due  to  the  functional 
dependence  of  the  charge-transfer  chemistry  (i.e.,  Eq.  (20))  and 
partially  due  to  compensating  temperature  effects.  At  low  flow 
rates,  the  MEA  temperatures  are  slightly  higher  owing  to  reduced 
convective  losses.  Because  of  high  activation  energies,  even  small 
temperature  variations  can  be  significant. 

Although  the  flow-rate  effects  on  polarization  may  be  weak,  it 
is  desirable  to  establish  the  minimum  flow  rates  needed  to  assure 
that  measured  results  are  independent  of  flow  rates.  Quantitative 
interpretation  of  experimental  results  depends  upon  quantitative 
knowledge  of  temperature  and  composition. 

7.  Performance  on  CH4-C02 

The  next  set  of  simulations  is  based  upon  a  fuel  mixture  of  75% 
CEI4  and  25%  C02  at  a  relatively  low  flow  rate  of  25  seem  (average 
velocity  of  0.16  m  s-1 ).  The  MEA  parameters  in  Table  1  continue  to 
be  used,  but  detailed  reforming  chemistry  is  now  included  [8].  The 
fuel  mixture  is  chosen  to  investigate  the  effects  of  internal  reform¬ 
ing.  In  this  case,  methane  reforming  proceeds  via  a  combination  of 
two  routes.  One  is  dry  reforming,  using  C02as  a  reforming  agent. 
Second  is  steam  reforming,  using  electrochemically  produced  EI20. 
Under  such  operating  conditions,  carbon  deposits  are  also  a  con¬ 
sideration.  Elowever,  the  present  study  does  not  directly  address 
carbon  deposition.  Because  reforming  is  an  endothermic  process, 
thermal  effects  within  the  MEA  are  more  complex  than  they  are 
with  H2fuel. 

Fig.  11  shows  contour  plots  of  species  mole  fractions  and  tem¬ 
perature  for  a  nominal  operating  condition  of  £cen  =  0.5  V.  With 
very  low  Reynolds  number  Re  (order  unity),  the  uniform  inlet  veloc¬ 
ity  relaxes  quickly  within  the  feed  tubes  to  the  laminar  Poiseuille 
parabolic  profile. 

Fig.  11  a  shows  contours  of  GT4(upper)  and  EI2  (lower)  mole 
fractions.  As  the  flow  approaches  the  MEA,  the  CH4  mole  fraction 
is  diluted  with  CO  and  H20.  Because  of  the  relatively  low  flow  rate, 
there  is  substantial  upstream  diffusion.  The  CEI4  mole  fraction  at  the 
anode  surface  is  approximately  half  of  the  specified  inlet  value.  The 
lower  plot  shows  significant  gas-phase  EI2  diffusion  upstream  into 
the  fuel  inlet  tube.  Because  of  the  low  flow  velocity  and  high  diffu- 
sivity  of  H2,  there  are  very  significant  variations  in  EI2  throughout 
the  anode  compartment. 
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Fig.  11.  Contour  plots  illustrating  the  solutions  for  the  cell  operating  with  75%  CH4, 
25%  CO2,  fuel  flow  rate  of  25  seem,  wall  temperature  of  800 °C,  and  Ece\\  =  0.5  V. 
The  composition  contours  are  mole  fractions.  Pathlines  are  superimposed  on  the 
contours  to  show  flow  patterns. 


Within  the  porous  anode,  CH4  is  catalytically  reformed  on  the 
Ni  catalyst  surfaces,  rapidly  producing  H2  and  CO.  At  the  interface 
with  the  dense  electrolyte  the  CH4  has  been  fully  converted.  The 
H2  reaches  a  peak  mole  fraction  of  around  30%  in  the  interior  of  the 
anode  structure.  Near  the  dense  electrolyte,  the  H2  is  consumed  by 
electrochemical  charge-transfer  reactions,  thus  significantly  reduc¬ 
ing  its  mole  fraction.  Relatively  high  CH4  mole  fractions  are  seen  in 
the  outer  radial  sections  of  the  MEA.  This  is  a  result  of  the  cathode 
being  applied  only  in  the  center  regions  of  the  disk,  eliminating 
electrochemistry  at  the  outer  radial  positions.  In  the  outer  regions 
the  electrochemistry  is  inactive  and  H20  is  not  formed,  which 
reduces  reforming  activity.  The  CH4  can  diffuse  radially  outward 
where  it  is  essentially  trapped.  In  a  similar  way  H2  can  diffuse  radi¬ 
ally  outward  within  the  anode  structure.  Because  H2  is  not  being 
produced  by  reforming  in  the  outer  radial  sections,  its  mole  fraction 
is  much  lower  than  in  the  central  regions  where  reforming  is  active. 

Fig.  11  b  shows  contours  of  C02(upper)  and  H20  (lower)  mole 
fractions.  The  electrochemically  produced  H20  is  maximum  within 
the  anode  near  the  dense-electrolyte  interface.  As  H20  is  trans¬ 
ported  toward  the  anode  compartment,  it  participates  in  catalytic 
reforming  chemistry  to  produce  H2.  Under  the  circumstances  here, 
the  H20  is  mostly  consumed  before  reaching  the  anode  compart¬ 
ment.  Significant  levels  of  C02  enter  with  the  fuel.  This  C02  is 
diluted  with  CO  and  H2  as  it  flows  toward  the  anode.  Within  the 
anode,  however,  C02  is  produced  catalytically,  primarily  via  water- 
gas-shift  chemistry.  That  is,  H20  reacts  with  CO  to  form  H2  and 
C02.  The  C02  mole  fraction  is  maximum  near  the  dense  electrolyte 
where  ample  H20  and  CO  are  present. 

Fig.  11  c  shows  mole  fraction  contours  of  CO  (upper)  and 
02( lower).  In  regions  where  charge-transfer  chemistry  is  active, 
the  H2  is  depleted  and  the  water-gas-shift  process  (globally,  CO  + 
H2 O  ^  C02  +  H2 )  rapidly  converts  CO  to  C02  and  H2 .  Thus  CO  levels 
remain  relatively  low  in  the  center  regions  of  the  MEA.  Interest¬ 
ingly,  the  CO  is  maximum  within  the  anode  structure  near  the 
outer  radius  of  the  cathode.  In  this  region,  CO  is  formed  primarily 
by  steam  reforming.  However,  because  of  reduced  electrochem¬ 
ical  activity  beyond  the  cathode  radius,  the  CO  is  not  converted 
into  C02  through  the  water-gas-shift  reaction  due  to  the  lower  H20 
concentration.  CO  also  diffuses  radially  outward,  where  it  is  essen¬ 
tially  trapped.  In  the  outer  regions  there  is  little  H20  available,  so 
water-gas  shift  is  relatively  inactive. 

On  the  cathode  side,  air  is  drawn  toward  the  cathode  where  it  is 
electrochemically  reduced  to  form  oxygen  ions  in  the  dense  elec¬ 
trolyte.  Thus,  the  02  mole  fraction  is  depleted  and  the  relative  N2 
fraction  increases.  In  this  system  there  is  not  a  forced  air  flow.  Nev¬ 
ertheless,  the  path  lines  show  induced  flow  as  air  is  drawn  toward 
the  cathode. 

Fig.  11  d  is  a  contour  plot  of  temperatures.  At  an  operating 
potential  of  £cell  =  0.5  V,  the  thermal  consequences  of  exothermic 
electrochemical  reactions  and  ohmic  heating  within  the  electrolyte 
are  nearly  balanced  by  the  endothermic  heterogeneous  reforming 
chemistry  within  the  anode.  Under  these  circumstances,  the  result 
is  nearly  isothermal  operation,  with  temperature  variations  of  only 
a  few  degrees.  Other  circumstances,  such  as  illustrated  in  Fig.  11  d, 
however,  cause  larger  temperature  excursions. 

Fig.  12  a  illustrates  the  methane  mole  fraction  as  a  function  of 
axial  position  on  the  button-cell  centerline.  The  position  x  =  0  is  at 
the  interface  between  the  porous  anode  and  the  fuel  gas  compart¬ 
ment.  At  the  relatively  low  fuel  flow  rate  of  25  seem,  it  is  evident 
that  there  is  significant  fuel  depletion,  especially  at  low  operating 
potential.  As  in  the  case  of  H2  fuel  (Fig.  7),  the  gradients  are  high¬ 
est  within  the  porous  electrode.  This  is  the  result  of  high  transport 
resistance  with  the  porous  media. 

Fig.  12  b  illustrates  H2mole  fraction  profiles.  Because  there  is 
no  H2  in  the  inlet  fuel  stream,  all  the  H2  is  produced  via  catalytic 
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Position  along  button-cell  axis  (mm) 

Fig.  12.  CH4and  H2  mole-fraction  profiles  along  the  button-cell  axis  within  the  fuel 
chamber  and  the  anode  structure.  In  all  cases,  the  inlet  fuel  composition  is  75%  CH4 
and  25%  C02  and  the  wall  temperature  is  800  °C. 

reforming  processes  within  the  anode.  The  H2  mole  fraction  is 
maximum  within  the  anode.  It  is  transported  toward  the  dense 
electrolyte  interface  where  it  is  consumed  via  charge-transfer 
reactions.  Hydrogen  is  also  transported  into  the  anode  gas  com¬ 
partment.  At  the  low  flow  rate  here  (25  seem),  the  H2  transport 
from  the  anode  structure  is  sufficient  to  maintain  significant  levels 
of  H2  within  the  anode  gas  compartment.  It  is  also  interesting  to 
note  that  at  high  operating  potentials  (low  current)  relatively  lit¬ 
tle  of  the  catalytically  formed  H2  is  consumed  electrochemically.  In 
this  case,  the  H2  is  formed  primarily  by  dry  reforming,  and  the  peak 
mole  fraction  is  near  the  dense  electrolyte  interface. 

Fig.  13  shows  radial  profiles  of  methane  mole  fraction  at  the 
interface  between  the  anode  and  the  fuel  chamber.  Significant 
radial  variations  are  evident.  It  is  interesting  to  note  that  the  CH4 
mole  fraction  increases  toward  the  outer  edges  of  the  anode.  Recall 
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Fig.  13.  Radial  profiles  of  CH4  mole  fraction  on  the  anode  top  surface.  In  all  cases, 
the  inlet  fuel  composition  is  75%  CH4  and  25%  C02.  The  wall  temperature  is  800  °C. 
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Fig.  14.  Temperature  profiles  along  the  button-cell  axis.  The  inlet  fuel  is  a  mixture 
of  75%  CH4  and  25%  C02  at  a  flow  rate  of  25  seem. 

that  in  this  cell,  the  outer  radius  of  the  cathode  layer  is  only  4  mm. 
Thus,  beyond  a  radius  of  4  mm,  there  is  no  electrochemical  steam 
production,  which  reduces  the  steam-reforming  rate.  Also,  some  of 
the  hydrogen  and  CO  produced  via  reforming  can  diffuse  radially 
into  the  outer  regions  of  the  porous  anode.  Because  in  this  region 
the  H2  is  not  consumed  electrochemically,  surface  chemistry  can 
regenerate  the  CH4.  The  simulations  here  suggest  that  a  higher  fuel 
flow  rate  is  desirable.  Increasing  the  flow  rate  will  reduce  radial 
variations  and  bring  the  fuel  composition  at  the  anode  face  closer 
to  the  inlet  composition. 

Fig.  14  shows  axial  temperature  profiles  along  the  button-cell 
axis.  In  all  cases  the  fuel  is  a  mixture  of  75%  CH4  and  25%  C02  at  25 
seem  and  the  furnace  walls  are  fixed  at  800  °C.  At  high  current  (low 
cell  potential)  the  MEA  temperature  exceeds  the  gas  temperatures 
on  both  the  fuel  and  air  sides.  This  is  because  the  ohmic  heating 
and  charge-transfer  exotherms  exceed  the  endothermic  effects  of 
on-anode  reforming.  At  low  current,  however,  the  situation  is  oppo¬ 
site.  Under  these  conditions  the  reforming  endotherm  dominates 
and  the  MEA  temperature  is  below  the  temperatures  in  the  gas 
compartments. 

8.  Flow  rates  and  cell  performance 

Understanding  fuel-depletion  effects  is  an  important  aspect  of 
understanding  SOFC  performance.  Operating  with  high  flow  rates 
is  not  in  conflict  with  investigating  fuel-depletion  effects.  As  is 
done  in  button-cell  experiments  by  Jiang  and  Virkar  [12],  the  fuel 
compositions  can  be  deliberately  premixed  to  simulate  various 
depletion  environments  that  occur  in  full  stacks  and  systems.  How¬ 
ever,  each  button-cell  experiment  is  run  under  flooded  conditions 
such  that  there  is  no  ambiguity  about  the  gas-compartment  com¬ 
position. 

All  other  things  being  equal,  fuel  flow  rates  should  be  sufficiently 
high  that  fuel  depletion  is  negligible  and  the  gas-phase  composi¬ 
tion  in  the  anode  compartment  is  essentially  that  of  the  feed  gas. 
In  principle  this  approach  simplifies  analysis  and  interpretation  of 
experimental  data.  However,  when  using  hydrocarbon  fuels  and 
Ni-based  anodes,  carbon  deposits  are  potentially  problematic.  As 
fuel  flow  rates  are  reduced,  the  hydrocarbon  fuel  is  diluted  with 
steam  and  C02,  and  hence  the  propensity  for  carbon  deposition  on 
Ni  surfaces  is  reduced. 

Although  using  low  fuel  flow  rates  may  complicate  analysis,  pre¬ 
venting  carbon  deposits  is  essential.  A  cell  that  fails  to  operate  due 
to  carbon  fouling  is  of  very  little  value.  With  low  fuel  flow  rates 
and  high  depletion  and  dilution,  it  should  be  recognized  that  the 
composition  of  the  fuel  seen  by  the  MEA  can  be  significantly  differ¬ 
ent  from  the  inlet  feed  composition.  To  be  quantitatively  useful, 
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models  that  are  used  to  interpret  the  data  must  accommodate 
the  gas-phase  composition  variations  within  the  anode  compart¬ 
ment. 

Fig.  11  shows  that  with  low  fuel  flow  rates  the  CH4  is  very  sig¬ 
nificantly  diluted  with  H20  and  C02  at  the  face  of  the  anode,  which 
contains  Ni.  Such  dilution  of  the  hydrocarbon  enables  deposit-free 
operation,  even  in  the  presence  of  Ni.  However,  it  may  also  be  noted 
that  there  are  alternatives  for  deposit-free  operation.  One  possibil¬ 
ity  is  to  use  electrode  materials,  such  as  doped  perovskites,  that  do 
not  promote  coke  formation  [18].  In  this  case,  high  fuel  flow  rate 
and  lower  depletion  is  practical. 

9.  Natural  convection  and  axisymmetry 

With  local  temperature  variations  and  a  horizontal  orientation 
of  the  cell,  it  is  reasonable  to  anticipate  that  natural  convection 
might  play  a  role  and  induce  three-dimensional  flow  fields.  How¬ 
ever,  for  the  conditions  studied  here,  the  results  show  that  natural 
convection  plays  a  negligible  role. 

The  potential  influence  of  natural  convection  behavior  can  be 
understood  in  the  context  of  a  dimensionless  groups.  Buoyant 
effects  are  typically  negligible  when 

<  1,  (25) 

Re2 

where  the  Grashof  and  Reynolds  numbers  are  defined  as 

Cr.Sa^,  R«.“  (26) 

vz  v 

In  these  definitions  g  is  the  gravitational  constant,  is  the 
volumetric  thermal-expansion  coefficient,  AT  is  a  characteristic 
temperature  difference,  U  is  a  characteristic  velocity,  L  is  a  char¬ 
acteristic  length  scale,  and  v  is  the  kinematic  viscosity.  For  this  cell 
configuration  the  characteristic  length  scale  for  the  Reynolds  num¬ 
ber  is  the  distance  between  the  inlet  tube  and  the  anode  face  (here 
3  mm)  and  the  characteristic  velocity  is  the  mean  inlet  velocity. 
Because  of  the  horizontal  orientation,  the  characteristic  length  scale 


Fig.  15.  CO(s)  surface  fractional  coverage  in  the  mid-anode  plane.  The  cell  is  operat¬ 
ing  with  a  fuel  flow  rate  of  25  seem,  nominal  temperature  of  800  °C,  and  Ecen  =  0.5  V. 
The  coverage  ranges  from  low  in  the  center  regions  (blue  contours)  to  high  at  the 
outer  radii  (red  contours).  (For  interpretation  of  the  references  to  colour  in  this  figure 
legend,  the  reader  is  referred  to  the  web  version  of  the  article.) 


for  the  Grashof  number  is  the  MEA  diameter.  Assuming  a  character¬ 
istic  temperature  difference  on  the  order  of  10  °C,  GrRe2  ^  1.  Thus, 
the  fact  that  the  three-dimensional  simulations  deliver  essentially 
axisymmetric  results  is  entirely  reasonable. 

As  an  example  that  illustrates  one  aspect  of  the  axisymme¬ 
try,  consider  surface  coverages  within  the  porous  anode  structure. 
Fig.  15  shows  contours  of  the  CO(s)  coverage  in  a  plane  through 
the  middle  of  the  anode  (£cell  =  0.5  V).  The  axisymmetry  is  evident 
from  the  bulls-eye  shape  of  the  contours.  Similar  results  are  found 
for  all  fluid  variables  at  all  flow  rates,  temperatures  and  operating 
potentials. 

An  axisymmetric  model  for  this  button  cell  would  be  quite  accu¬ 
rate  (i.e.,  three-dimensional  effects  are  found  to  be  negligible).  Of 
course,  larger  dimensions,  lower  flow  rates,  or  higher  temperatures, 
can  lead  to  circumstances  in  which  buoyant  effects  could  become 
important. 

10.  Summary  and  conclusions 

Button-cell  experiments  are  used  widely  in  the  development 
of  new  SOFC  materials  and  MEA  architectures.  Data  usually  con¬ 
sist  of  polarization  characteristics  and  sometimes  electrochemical 
impedance  spectra.  The  experimental  conditions  are  usually  rep¬ 
resented  in  terms  of  inlet  fuel  composition  and  flow  rate,  cell 
temperature,  and  operating  potential. 

In  analyzing  and  interpreting  button-cell  experiments,  two 
important  assumptions  are  typically  applied.  First,  it  is  assumed 
that  the  flow  rates  are  sufficiently  high  that  the  species  composi¬ 
tion  at  the  interface  between  the  gas  compartment  and  the  anode 
structure  is  the  same  as  inlet  fuel  mixture.  In  other  words,  the  fuel 
is  minimally  depleted  and  diluted.  Second,  the  gases  and  the  MEA 
are  assumed  to  be  isothermal.  If  either  of  these  assumptions  is  not 
valid,  then  the  data  interpretation  can  be  flawed. 

In  this  paper  a  new  three-dimensional  computational  fluid 
dynamics  model  is  used  to  characterize  physical  and  chemical 
behavior  in  SOFC  button  cells.  The  primary  purpose  is  to  investi¬ 
gate  how  the  assumptions  noted  above  depend  upon  cell  operating 
conditions.  The  model  includes  coupled  fluid  flow,  heat  transfer, 
and  multi-component  species  transport.  It  also  includes  porous- 
media  transport,  elementary  catalytic  chemistry  within  electrode 
structures,  and  electrochemical  charge  transfer. 

The  model  is  implemented  in  object-orientated  C++.  It  uses 
geometrically  flexible  unstructured  grids,  which  incorporate  fluid 
regions  and  solid  materials.  The  open-source  Cantera  software  is 
used  to  evaluate  thermodynamic  and  transport  properties,  as  well 
as  reaction  chemistry. 

Although  the  model  can  be  applied  to  general  geometries,  it  is 
used  here  to  simulate  a  particular  button  cell  that  has  been  exper¬ 
imentally  measured.  Results  show  that  if  inlet  flow  rates  are  too 
low,  there  can  be  significant  compositional  variations  in  the  gas 
compartments.  When  this  is  the  case,  it  is  not  appropriate  to  ana¬ 
lyze  data  using  the  assumption  that  the  MEA  is  exposed  to  the  fuel 
composition  that  is  supplied  through  the  feed  lines.  However,  it  is 
possible  to  increase  flow  rates  to  the  point  that  gas  composition 
is  minimally  depleted  and  diluted.  Models,  such  as  the  one  in  this 
paper,  provide  the  quantitative  means  to  establish  the  appropriate 
operating  conditions. 

Because  the  button-cell  in  this  study  is  oriented  horizontally, 
there  is  the  possibility  for  buoyant  convection  to  cause  three- 
dimensional  flow  patterns.  Such  behavior  could  further  complicate 
data  interpretation.  However,  because  of  relatively  small  dimen¬ 
sions  and  temperature  variations,  the  results  show  that  buoyant 
convection  is  easily  overcome.  Thus,  even  under  circumstances 
with  low  flow  rates  high  fuel  depletion,  assuming  axisymmetric 
flow  is  valid. 


G.M.  Goldin  et  al.  /  Journal  of  Power  Sources  187  (2009)  123-135 


135 


Acknowledgements 

This  work  at  the  Colorado  School  of  Mines  and  Ansys/Fluent  was 
supported  by  the  Office  of  Naval  Research  through  a  Research  Tools 
Consortium  grant  number  N00014-05-1-0339.  The  effort  at  North¬ 
western  University  was  supported  by  the  Department  of  Energy 
National  Energy  Technology  Laboratory  under  Award  Number  DE- 
FC26-05NT42625.  We  gratefully  acknowledge  the  collaboration 
with  Prof.  David  Goodwin  (Caltech)  for  his  assistance  in  integrating 
Cantera  into  these  simulations. 

References 

[  1  ]  Z.  Zhan,  Y.  Lin,  M.  Pillai,  I.  Kim,  S.A.  Barnett,  J.  Power  Sources  161  (2006)  460-465. 

[2]  Y.  Shi,  N.  Cai,  C.  Li,  J.  Power  Sources  164  (2007)  639-648. 

[3]  V.M.  Janardhanan,  O.  Deutschmann,  J.  Power  Sources  162  (2006)  1192-1202. 

[4]  R.  Suwanwarangkul,  E.  Croiset,  E.  Entchev,  S.  Charojrochkul,  M.D.  Pritzker,  M.W. 
Fowler,  P.L.  Douglas,  S.  Chewathanakup,  H.  Mahaudom,  J.  Power  Sources  161 
(2006)308-322. 

[5]  Y.  Shi,  N.  Cai,  C.  Lia,  C.  Bao,  E.  Croiset,  J.  Qian,  Q.  Hu,  S.  Wang,  J.  Power  Sources 
172 (2007) 235-245. 


[6]  Y.  Shi,  N.  Cai,  C.  Lia,  C.  Bao,  E.  Croiset,  J.  Qian,  Q.  Hu,  S.  Wang,  J.  Power  Sources 
172(2007)246-252. 

[7]  D.S.  Monder,  K.  Nandakumar,  K.T.  Chuang,  J.  Power  Sources  162  (2006)  400- 
414. 

[8]  H.  Zhu,  R.J.  Kee,  V.M.  Janardhanan,  O.  Deutschmann,  D.G.  Goodwin,  J.  Elec- 
trochem.  Soc.  152  (2005)  A2427-A2440. 

[9]  H.  Zhu,  R.J.  Kee,  J.  Electrochem.  Soc.  155  (2008)  B715-B729. 

[10]  R.J.  Kee,  M.E.  Coltrin,  P.  Glarborg,  Chemically  Reacting  Flow:  Theory  and  Practice, 
John  Wiley,  New  York,  2003. 

[11  ]  H.  Zhu,  R.J.  Kee,  ECS  Trans.  7  (2007)  1869-1878. 

[12]  Y.  Jiang,  A.V.  Virkar,  J.  Electrochem.  Soc.  150  (2003)  A942-A951. 

[13]  A.M.  Sukeshini,  B.  Habibzadeh,  B.P.  Becker,  C.A.  Stoltz,  B.W.  Eichhorn,  G.S.  Jack- 
son,  J.  Electrochem.  Soc.  153  (2006)  A705-A715. 

[14]  S.  Mathur,  J.  Murthy,  Numer.  HeatTransf.  31  (1997)  195-215. 

[15]  S.V.  Patankar,  Numerical  Heat  Transfer  and  Fluid  Flow,  McGraw-Hill,  New  York, 
1980. 

[16]  K.  Stueben,  in:  C.W.  Oosterlee,  U.  Trottenberg,  A.  Schuller  (Eds.),  Multigrid, 
Academic  Press,  New  York,  2001,  pp.  413-532. 

[17]  D.G.  Goodwin,  in:  M.  Allendorf,  F.  Maury,  F.  Teyssandier  (Eds),  Chemical  Vapor 
Deposition  XVI  and  EUROCVD  14,  vol.  PV  2003-08,  pp.  155-162.  Electrochem¬ 
ical  Society,  2003  (see  also  http://www.cantera.org). 

[18]  M.R.  Pillai,  Y.  Jiang,  N.  Mansourian,  I.  Kim,  D.M.  Bierschenk,  H.  Zhu,  R.J.  Kee,  S.A. 
Barnett,  Electrochem.  Solid-State  Lett.  11  (2008)  B174-B177. 


